%% RK4
clc; clear; close all;

% 参数
alpha = 2.4;
A = [-alpha, 0, 0;
     0, 0, 0;
     0, 0, 0];
B = [1, -4, -3.5;
     0, 1, 2;
    -1, -4, 1.5];

% 只跑这一个q10
q10_list = 0.0623;  

tspan = [0 2000];
dt = 0.01;
N = floor((tspan(2) - tspan(1)) / dt);

f = @(t,x) mCNN_single(t, x, A, B);

results = cell(length(q10_list),1);

parfor idx = 1:length(q10_list)
    q10 = q10_list(idx);
    x0 = [1e-6; 0; 0; q10];
    
    X = zeros(length(x0), N+1);
    X(:,1) = x0;
    
    for i = 1:N
        X(:,i+1) = RK4(f, X(:,i), dt);
    end
    
    results{idx}.X = X;
    results{idx}.q10 = q10;
end

% 画图
T = linspace(tspan(1), tspan(2), N+1);
sel = T >= 1000;
X_sel = results{1}.X(:, sel);

figure;
plot(X_sel(1,:), X_sel(3,:), 'LineWidth', 2);
xlabel('x_1');
ylabel('x_3');
title(sprintf('相空间轨迹 (q_{10}=%.3f)', results{1}.q10));
% xlim([-3, 3]);
% ylim([-1.5, 1.5]);
grid on;

% 微分方程函数
function dxdt = mCNN_single(~, x, A, B)
    X = x(1:3);
    q1 = x(4);

    M = exp(-q1);
    Y = M .* X;

    dX = -X + A*Y + B*X;
    dq1 = X(1);

    dxdt = [dX; dq1];
end


%% ODE45
% clc; clear; close all;
% 
% % 参数
% alpha = 2.4;
% A = [-alpha, 0, 0;
%      0, 0, 0;
%      0, 0, 0];
% B = [1, -4, -3.5;
%      0, 1, 2;
%     -1, -4, 1.5];
% 
% % 初始条件
% x0 = [1e-6; 0; 0];  % 状态变量初始
% q10 = 0.9;            % 内部变量初始
% 
% % 时间设置
% tspan = [0 2000];
% dt = 0.01;
% t_eval = 0:dt:tspan(2);
% 
% % 定义微分方程
% odefun = @(t,x) mCNN_single(t, x, A, B);
% 
% % 初始条件扩展，包含 q1
% x_init = [x0; q10];
% 
% % 数值积分
% [T, X] = ode45(odefun, t_eval, x_init);
% 
% % 取稳态段绘图
% sel = T >= 50;
% X_sel = X(sel,:);
% 
% figure;
% plot(X_sel(:,1), X_sel(:,3), 'LineWidth', 2);
% xlabel('x_1');
% ylabel('x_3');
% title(sprintf('相空间轨迹 (q_{10}=%.3f)', q10));
% % xlim([-3  , 3  ]);
% % ylim([-1.5, 1.5]);
% grid on;
% 
% function dxdt = mCNN_single(~, x, A, B)
%     X = x(1:3);
%     q1 = x(4);
% 
%     M = exp(-q1);
%     Y = M .* X;
% 
%     dX = -X + A*Y + B*X;
%     dq1 = X(1);
% 
%     dxdt = [dX; dq1];
% end